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I— ABSTRACT 

Context. Previous work has shown that the tidal interaction between a binary system and a circumbinary disc leads to the formation of a large 
Q inner cavity in the disc. Subsequent formation and inward migration of a low mass planet causes it to become trapped at the cavity edge, where 
^ ' it orbits until further mass growth or disc dispersal. The question of how systems of multiple planets in circumbinary discs evolve has not yet 
C/^ , been addressed. 

■ Aims. We present the results of hydrodynamic simulations of multiple low mass planets embedded in a circumbinary disc. The aim is to 
examine their long term evolution as they approach and become trapped at the edge of the tidally truncated inner cavity. 

Methods. A grid-based hydrodynamics code was used to compute simulations of 2D circumbinary disc models with embedded planets. The 
^ 3D evolution of the planet orbits was computed, and inclination damping due to the disc was calculated using prescribed forces. We present a 
^—H _ suite of simulations which study the evolution of pairs of planets migrating in the disc. We also present the results of hydrodynamic simulations 

■ of five-planet systems, and study their long term evolution after disc dispersal using a N-body code. 

Results. For the two-planet simulations we assume that the innermost planet has migrated to the edge of the inner cavity and remains trapped 
there, and study the subsequent evolution of the system as the outermost planet migrates inward. We find that the outcomes largely depend 
fSj ' on the mass ratio q = mi/rrio, where (m^) is the mass of the innermost (outermost) planet. For q < I, planets usually undergo dynamical 
^— H " scattering or orbital exchange. For values of ^ > 1 the systems reach equilibrium configurations in which the planets are locked into mean 
motion resonances, and remain trapped at the edge of the inner cavity without further migration. Most simulations of five-planet systems we 
performed resulted in collisions and scattering events, such that only a single planet remained in orbit about the binary. In one case however, a 
^ multiplanet resonant system was found to be dynamically stable over long time scales, suggesting that such systems may be observed in planet 
searches focussed on close binary systems. 
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1. Introduction 

At the time of writing approximately 250 extrasolar planets 
have been discovered, of which about 30 are members of binary 
or multiple stellar systems (Eggenberger et al. 2004; Mugrauer 
et al. 2007). Most of the planets found in binaries orbit around 
one of the stellar component in so-called S-type orbits, and 
the majority of binaries harbouring planets have orbital separa- 
tions Gh > 100 AU. There are, however, exceptions to these 
cases. The short period binary systems Gliese 86, y Cephei 
and HD41004 have at ~ 20 AU, and contain planets orbit- 
ing at 1 - 2 AU from the central star (Eggenberger et al. 2004; 
Mugrauer & Neuhauser 2005). Although there are no known 
planets which orbit around both stellar companions in a bi- 
nary system consisting of main sequence stars (i.e. circumbi- 
nary planets on so-called P-type orbits), there are two sys- 
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tems which indicate that the formation of circumbinary plan- 
ets is feasible. The first is the circumbinary planet with mass 
nip = 2.5 M/ and orbital radius 23 AU observed in the radio 
pulsar PSR 1620-26 (Sigurdsson et al. 2003). The second is the 
nip = 2.44 Mj planet orbiting about the system composed of 
the star HD202206 and its 17.4 M/ brown dwarf companion 
(Udry et al. 2002). The lack of observed circumbinary planets 
is probably due to the fact that short period binaries are usually 
rejected from observational surveys. 

The observation of planets in binary systems is consis- 
tent with detections of circumstellar discs in binary systems. 
Several circumbinary discs have been detected around spectro- 
scopic binaries such as DQ Tau, AK Sco, and GW Ori. The 
circumbinary disc around GG Tau has been directly imaged 
(Dutrey et al. 1994), revealing the presence of a tidally trun- 
cated inner cavity generated by the central binary. The exis- 
tence of these circumbinary discs opens up the possibility of 
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circumbinary planets forming. Combining the observations of 
circumbinary discs, with the fact that ~ 50% of solar-type stars 
are members of binaries (Duquennoy & Mayor 1991), suggests 
that circumbinary planets are probably common. 

To date there has been a relatively modest amount of the- 
oretical work examining planet formation in binary systems. 
Results from previous studies suggest that planetesimal ac- 
cretion should be possible in regions of both circumstellar 
(Marzari & Scholl 2000; Thebault et al. 2006) and circumbi- 
nary discs (Moriwaki & Nakagawa 2004; Scholl et al. 2007). 
Quintana & Lissauer (2006) simulated the late stages of ter- 
restial planet formation in circumbinary discs. They found that 
planetary systems similar to those around single stars can be 
formed around binaries, provided the ratio of the binary apoc- 
entre distance to planetary orbit is < 0.2. In general binaries 
with larger maximum separations lead to planetary systems 
with fewer planets. 

The evolution of a low mass planetary core embedded in 
a circumbinary disc was investigated recently by Pierens & 
Nelson (2007) (hereafter referred to as Paper I). This work ex- 
amined the migration and long term orbital evolution of planets 
with masses of =5, 10 and 20 under the action of disc 
torques. It was found that the inward drift of a planet undergo- 
ing type I migration is stopped at the edge of the cavity formed 
by the binary. This halting of migration is due to positive coro- 
tation torques operating which can counterbalance negative 
Lindblad torques. Such an effect is known to be at work in ac- 
cretion disc regions where there is a strong positive gradient of 
the surface density (Masset et al. 2006). Interestingly, Pierens 
& Nelson (2007) showed that the stopping of migration in cir- 
cumbinary discs occurs in a region of long-term dynamical sta- 
bility, suggesting that such planets may be able to survive there 
over long times, or at least remain in the disc for long enough 
to form a gas giant planet. The evolution of giant planets in 
circumbinary discs was considered by Nelson (2003). 

In this paper, we extend the model presented in Paper I by 
considering the evolution of multiple planets embedded in a cir- 
cumbinary disc. Here, we wish to examine how multiple plan- 
ets interact with each other if they form at large distance from 
the binary and successively migrate toward the cavity edge. In 
particular, we want to look at whether or not the trapping of 
a planet at the cavity edge, and the subsequent migration of 
additional planets to its vicinity lead to growth of the planet 
through collisions, the formation of mean motion resonances, 
or destablisation of the system through gravitational scattering. 

To address these issues, we first consider a system which 
consists of a pair of planets with masses of = 5, 10 and 20 
M®. We assume that one planet is trapped at the edge of the cav- 
ity while the outermost planet migrates in from larger radius. 
The results of the simulations show that the final outcome of 
such a system generally depends on the mass ratio q = mi/rrio 
(where m/ is the mass of the inner planet and is the mass of 
the outer planet). Interestingly, we find that systems with q> 1 
can reach a steady state such that the planets are locked into 
resonance and remain trapped at the cavity edge. Most of the 
systems with q < I, however, are unstable and lead to events 
such as scattering or dynamical exchange. 



We performed a second set of simulations consisting of 
five-planet systems embedded in a circumbinary disc. Of three 
simulations performed, two resulted in a single planet orbiting 
around the binary because of collisions and scattering events. 
The remaining simulation resulted in a three-planet system re- 
maining, with all planets in mutual mean motion resonances. 
This configuration was found to be stable over long time scales. 

This paper is organized as follows. In section 2 we describe 
the physical model and the numerical method. In section 3 we 
describe the results of simulations aimed at studying the evo- 
lution of pairs of planets embedded in a circumbinary disc. We 
then present in section 4 the simulations of five-planet systems 
embedded in a circumbinary disc. We finally discuss our results 
and present our conclusions in section 4. 

2. Physical model and numerical setup 

2.1. Disc and planet evolution 

As in Paper I, we adopt a 2D disc model for which we assume 
no vertical motion. The equations governing the disc evolution 
are described in detail in Paper I and therefore will not be dis- 
cussed here. 

In Paper I, the planet orbit and the disc midplane were assumed 
to be coplanar. However, the simulations presented here exam- 
ine the evolution of multiple planets which can strongly interact 
with each other as their orbits converge, leading eventually to 
close encounters. During such an event, a planet may receive a 
significant component of acceleration in the vertical direction, 
reducing thereby the interaction with both the disc and other 
planets. As a consequence, we decided here to use a model 
in which planets can evolve in the z direction as well. With 
respect to coplanar orbits, this will also reduce the collision 
rate between planets, increasing thereby the time during which 
planets can strongly interact. 

In the work presented here, each planet can interact with every 
other one and with the disc. The latter interaction is expected 
to lead not only to orbital migration but also to eccentricity and 
inclination damping. The gravitational potential of the disc is 
calculating using the following expression: 



[ 



bp) + 4 + 6^ 



(1) 



where S is the disc surface density and where the integral is 
performed over the disc surface, r^, 0^ and Zp are respectively 
the radial, azimuthal and vertical coordinates of the p^^ planet. 
6 is a smoothing parameter which is set to 6 = 0.6 H, where H is 
the local disc scale height. Under the action of this gravitational 
potential, each planet undergoes both orbital migration and ec- 
centricity damping. However, because of the 2D disc model we 
use here, bending waves cannot be launched in the disc, and so 
there is no disc induced damping of inclination. To model the 
latter we follow Tanaka & Ward (2004) and mimic the eff'ect of 
bending waves by applying to each planet a vertical force 
given by: 
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where Cg is the sound speed and where Q and Hp are respec- 
tively the Keplerian angular velocity and the disc surface den- 
sity at the position of the planet. and are dimensionless 
coefficients which are set to = -1.088 and = -0.871 
(Tanaka & Ward 2004), and is a free parameter which is 
chosen such that the inclination damping ti timescale obtained 
in the simulations is approximatively equal to the eccentric- 
ity damping timescale tg. Test simulations show that choosing 
13 = 0.33 give similar values for U and tg. We adopt therefore 
this value for this work. Here, it is worthwhile to notice that ac- 
cording to the linear theory (Tanaka & Ward 2004), a small 
value of the planet inclination (ip «c H/R) is not expected 
to affect the migration rate of the planet. For large values of 
ip however, migration rates may be moderately slowed down 
(Cresswell et al. 2007) because the interaction with the disc 
is reduced. Such an effect is accounted for in an approximate 
manner by Eq. [T] 

2.2. Numerical setup 
2.2.1. Numerical method 

The simulations presented here were performed using the hy- 
drodynamic code GENESIS. This code employs a second- 
order- accurate method that computes advection using the 
monotonic transport algorithm (Van Leer 1977). Details about 
GENESIS are given in Paper I. All the runs use Nr = 256 radial 
grid cells uniformly distributed between r/„ = 0.5 and rout = 6 
and A/^^ = 380 azimuthal grid cells. 

The evolution of the planets plus binary system is performed 
using a 5th-order Runge-Kutta scheme (Press et al. 1992). In 
spite of the accuracy of this integrator, experiments have shown 
that to ensure energy conservation during close encounters, the 
time step size used to make the system evolve should be 
smaller than the hydrodynamical time step Ath based on the 
GEL criterion (Stone & Norman 1992). Eollowing Gresswell 
& Nelson (2006), we ensure good energy conservation and ac- 
curacy by setting the time step size to At=mm(Ath, Atn\,Atni), 

where: 

In 
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In the previous expressions, Vpp' is the distance between the 
planets p and p' and r^^ is the distance between the planet p 
and star s. We note that throughout our simulations, the time 
step size is ~ 1 /800 the binary orbital period. 
As in Paper I, we adopt computational units in which the to- 
tal mass of the binary is M., = 1 , the gravitational constant is 
G = 1 , and the radius r = 2 in the computational domain cor- 
responds to 5 AU. The unit of time is Q"^ = ^GM^/a^, where 
at = OA is the initial semimajor axis of the binary. This corre- 
sponds to an initial separation between the two stars of ~ 1 AU. 
In the simulations presented here, close encounters between 
two planets can result in a physical collision. Here, this is 



planets p and p' is less than (3mp/4npy'^ -h (3mp/47ipy 
where p is the mass density which we assume to be the same for 
each planet and equal to p = 3 g.cnT^. If a collision is found 
to occur between the planets p and p\ these are assumed to 
merge and are subsequently substituted by a single body with 
mass nip -\- nip. The position and velocity of the latter are set to 
the position and velocity of the centre of mass of the planets p 
and p'. 



2.2.2. Initial conditions 

As in Paper I, the disc aspect ratio H/R is assumed to be 
constant and equal to H/R = 0.05. We use also the "alpha" 
prescription of Shakura & Sunyaev (1973) to model the disc 
anomalous kinematic viscosity v = aCsH, where Cs is the 
isothermal sound speed and where a = 10~^. The reason for 
choosing such a low a value is discussed in detail in Paper I, 
but is essentially because a larger value causes rapid evolution 
of the binary orbit that would prohibit the long simulations we 
present here. 

In Paper I, we showed that both the disc and binary evolve to- 
ward a near- steady state as they interact with each other. Erom 
the time this equilibrium configuration is reached, the apsidal 
lines of the disc and binary are aligned. Also, the disc struc- 
ture and the orbital elements of the binary remain essentially 
constant. Eor example, we find that the eccentricity of a binary 
with mass ratio = 0.1 and initial separation at = 0.4 satu- 
rates at a value of ~ 0.08. The simulations presented in paper 
I of one planet interacting with a circumbinary disc were per- 
formed using this quasi-equilibrium state as initial conditions 
for the disc and binary. Depending on the model we consider, 
we adopt here a similar approach when setting up our initial 
conditions: 

i) In simulations of pairs of planets embedded in a circumbi- 
nary disc, we restart the runs presented in Paper I once the 
planet is trapped at the cavity edge but with a second planet 
evolving on a circular orbit with ap = 2.5 and ip = 0.5°. 
The latter is allowed to interact with the disc whose mass is 
Md -0.01 Mi, and with the other planet and binary. 

ii) In simulations that evolve five-planet systems in a circumbi- 
nary disc, we embed the planets in the disc once the latter and 
the binary have reached a stationary state, as described in Paper 
I. We set the innermost planet at = 1.8 and then calculate 
the initial location of the others by asssuming that two adjacent 
bodies p and p' are separated by ~ 5 RmH, where RmH is the 
mutual Hill radius defined by: 



R. 



mH 



3M^ 



1/3 



^Qp j 



(5) 



Each body is assumed to initially evolve on a circular orbit with 
ip = 0.5°. Note that the initial separation between planets we 
adopt is greater that the critical value of ~ 3.46 R^h below 
which rapid instability occurs for two planets on initially circu- 
lar orbits (Gladman 1993). 
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Table 1. The first column gives the run label, the second col- 
umn gives the mass m/ of the inner planet and the third column 
gives the mass of the outer planet. The fourth column gives 
the ratio q = mi/mo. 



Run jfij (Me) (M^) q 



Runl 


10 


10 


1 


Run2 


10 


5 


2 


Run3 


20 


10 


2 


Run4 


20 


5 


4 


Run5 


5 


10 


0.5 


Run6 


10 


20 


0.5 


Run? 


5 


20 


0.25 



3. Evolution of pairs of planets embedded in 
circumbinary discs 

In paper I we considered planets with masses mp= 5, 10 and 
20 M®. Therefore, we adopt here the same values for the mass 
Mi of the innermost planet which is assumed to be trapped at 
the cavity edge of the disc. For each value of m/, we have per- 
formed two or three runs for which ruo was varied between 
5 < < 20 M®. Table [T] gives the values of m/, and 
q = Mi I Mo for each run. An interesting feature of the results 
of these simulations is that varying the value of q can lead to 
diff'erent outcomes. Below, we discuss in detail the diff'erent 
modes of evolution obtained in the simulations, and how they 
change depending on whether or not q>\. 

3.1. Models with q = I 

The orbital evolution of a pair of planets with masses of ntp = 
10 M© is shown in the left panel of Figure [H At the beginning 
of the simulation, corotation torques cause the inner planet to 
remain trapped at the edge of the disc cavity (see Paper I) lo- 
cated at r ~ 1.2. The outer planet, which initially evolves on a 
circular orbit with Gq = 2.5, undergoes the usual type I migra- 
tion and drifts inward. As it migrates, its eccentricity eo slowly 
increases due to the influence of the binary (see bottom left 
panel of Figure [T]). Also, the distance between the two bod- 
ies becomes smaller and smaller which can potentially leads 
to the formation of mean motion resonances (Papaloizou & 
Szuszkiewicz 2005). Here, we find that the two planets become 
captured into the 4:3 resonance at ^ ~ 1.9 x 10^ binary or- 
bits. The top right panel of Figure [T] display s the evolution of 
both the apsidal angle Aoj = oji - ojo and the resonant angle 
i// = 4Ao - 3Ai - 3aji, where Ai (Ao) and oJi (oJo) are respec- 
tively the mean longitude and longitude of pericentre of the in- 
ner (outer) planet. Once the resonance is established, the libra- 
tion amplitude of slightly increases with time until ^ ~ 4x 10^ 
binary orbits, and then remains almost unchanged, suggesting 
that the planets are stably locked into the resonance. From this 
time onward, the system is close to an equilibrium state with 
both planets having constant semimajor axes and eccentrici- 
ties. At the end of the simulation, the ratio of semimajor axes 
is Gi/ao ~ 0.85 and the ratio of eccentricities is ei/eg ~ 0.6. 
Such a configuration can be achieved because the torques ex- 
erted by the disc on each planet act in an opposite way and can 



eventually counterbalance each other. From the time the res- 
onance is established the negative torques exerted by the disc 
on the outer planet make the two planets migrate inward to- 
gether. However, as both planets migrate, the innermost planet 
experiences stronger positive corotation torques which tend to 
push the pair of planets outward. The bottom right panel of Fig. 
\T\ shows the evolution of the torques exerted on each planet as 
well as the eff'ective torques acting on the whole system. We see 
that as the evolution proceeds, the torques exerted on the inner 
planet are able to exactly counterbalance the ones exerted on 
the outer body, which leads consequently to a zero net torque 
acting on the system. This happens from r ~ 2.1 x 10^, thereby 
stopping the joined migration of the planets. 



3.2. Models with q > 1 

Stable resonant locking was also found in some of the calcu- 
lations with q > I. Fig. [2] shows the results for Run2 in which 
planets have masses of m/ = 10 and rUo = 5 M©. With respect 
to Runl, the outer planet migrates more slowly since its mass 
is smaller. However, the mode of evolution found in Run2 is 
very similar to the one obtained in Runl, leading ultimately 
to a stable configuration with the two planets trapped in the 
4:3 resonance from ^ ~ 5 x 10"^. At earlier times, the evolu- 
tion of Co shows some peaks at ^ ~ 1.9 x 10"^ and r ~ 3 x 10"^ 
which coincide with the planets being temporary captured in 
the 2:1 and 3:2 resonances. At the end of the simulation, Co is 
still slightly increasing whereas Ci is slightly decreasing, which 
indicates that the equilibrium configuration is not fully estab- 
lished. Nonetheless, comparing Figs. [T] and |2l we can see that 
the libration amplitude of the resonant angle is much smaller in 
Run2 than in Runl, suggesting that the 4:3 resonance is more 
stable in this case. 

In models with m/ = 20 M©, the simulations resulted in dif- 
ferent modes of evolution, depending on the mass of the outer 
planet m^. Fig. [3] and Fig. |4l display the results of calculations 
with nto = 10 M© and rrio = 5 M© respectively. Comparing 
these two figures, we can see that the final state of the system 
is quite similar in both cases, with planets evolving on fixed 
orbits with ai ~ 1.2 and Go ~ 1.6. In the run with = 10 M©, 
a 3:2 resonance forms at ^ ~ 1.4 x 10^. In the simulation with 
nto = 5 M© however, there is no evidence that the planets are 
in mean motion resonance, even though the system is close 
to the 3:2 commensurability. In this case, examination of the 
torques exerted by the disc (see upper panel of Fig. [5]) reveals 
that the migration of the system stalls because the torques act- 
ing on both planets cancel. Here, such an eff'ect arises because 
the mass of the inner planet is high enough to make the on- 
set of non-linear eff'ects possible. These can significantly alter 
the surface density profile and widen the size of the inner cav- 
ity. Indeed, we find that the edge of the inner cavity is located 
at r ~ 1.2 in simulations with rui = 10 M© and m/ = 5 M© 
whereas the bottom panel of Fig. [5] shows that it is located at 
r ~ 1.5 in Run4. Consequently, the evolution of the system in 
Run4 is such that the migration of the outermost planet is halted 
at the edge of the cavity formed by the binary plus inner planet 
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Runi A:S r-esonance 




^>t^u* 2x10."* 5x1 4xid^ 2^\o* jmld^ 4x10:^ 

Time (tjlnqry ^rbUs-) Time ^tilnqry q-rblts) 



Fig. 1. Left: evolution of the semimajor axes (top panel) and eccentricities (bottom panel) of planets with masses of = 10 
Me. Right: the top panel shows the evolution of both the resonant angle i// = 4Ao - 3Ai - Soji associated with the 4:3 resonance 
(black) and the apsidal angle Aoj = oji - cOq (green). The bottom panel shows the evolution of the torques exerted by the disc on 
the innermost planet (red line) and on the outermost one (yellow line). The sum (red+yellow) is displayed with the black line and 
corresponds to the torques exerted on the whole system. 



system, therefore preventing capture in the 3:2 resonance (or in 
resonances of higher degree such as 4:3). 

3.3. Models with q < 1 

For simulations with q < 1, the evolution of the system dif- 
fered significantly from that just described in all but one case 
(Run6). Fig. [6] shows the results of a simulation (Run5) with 
q = 0.5 in which planets have masses of m/ = 5 M© and 
rrio = 10 M©. Here, the outer planet passes through the 4:3 
resonance at r ~ 1.8 x 10"^ and then slips into the 6:5 reso- 
nance with the inner planet at ^ ~ 2 x 10"^. The passage through 
these resonances is clearly accompanied by an increase of the 
inner planet eccentricity At ^ ~ 2.1 x 10"^, the inner planet 
undergoes a close encounter with the outer one as a result of 
resonant trapping. This subsequently leads to the scattering of 
the inner planet further out in the disc while the outer one is 
pushed inward by virtue of conservation of angular momen- 
tum. Interestingly, such an orbital exchange leads to a config- 
uration of the system similar to that of models with q > 1 . In 
good agreement with what we described in Section [T2l we find 
that the final state of the system is indeed an equilibrium con- 
figuration with the two planets locked into the 6:5 resonance. 
At the end of the simulation, the 5 M© and 10 M© planets are 
respectively located at a/ ~ 1.35 and Go ~ 1.2. 
Although q has the same value in Run6 as it is in Run5, we find 
a diff'erent mode of evolution as shown in Fig. [71 In Run6 the 
planets have masses jtii — 10 Af© and jtIq — 20 Af©. Relative 
to Run5 the disc induced eccentricity damping acting on the 
innermost planet is stronger, thereby preventing eccentricities 
from reaching large values, and consequently preventing plan- 
ets from undergoing close encounters. Fig. [7] shows that the 
mode of evolution of the system in Run6 is almost similar to 



that of models with q > I, since the planets are in resonance 
(the 6:5 resonance in this case) and do not migrate. At the end 
of the evolution, the innermost and outermost planets are re- 
spectively located at a/ ~ 1.1 and Gg ~ 1.3 and the ratio of 
eccentricities is ei/eg ~ 1.3. 

Fig.[8]shows the results of a calculation (Run?) with nti = 5 M© 
and nto = 20 M©, corresponding to a model with q = 0.25. 
Here, the 20 M© planet enters the 4:3 resonance with the 5 M© 
body at r ~ 9 x 10^, which drives the eccentricity of the latter 
upward. This occurs until the inner planet has a close encounter 
with the binary, resulting in the planet being completely ejected 
from the system. At later times, the evolution is close to that de- 
scribed in Paper I, with the 20 M© planet migrating until it is 
trapped at the edge of the cavity formed by the binary. 

4. Evolution of five-planet systems embedded in 
circumbinary discs 

We now turn to the question of how a swarm of planets em- 
bedded in a circumbinary disc evolves. To address this issue, 
we have considered diff'erent models in which five planets hav- 
ing masses of 5, 7.5, 10, 12.5 and 15 M© interact with each 
other. Three simulations have been performed in which we var- 
ied the initial configuration of the system. In one run (hereafter 
Model 1), the initial mass distribution of planets decreases as 
one moves out in the disc, whereas it is chosen to be random in 
Model2. In Model3 the initial mass distribution increases as a 
function of increasing orbital radius. 

4.1. Modell 

In all the simulations we have performed, the innermost planet 
initially evolves on a circular orbit with ai = 1.8. For this 
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Fig. 2. This figure shows the evolution of the semimajor axes 
and eccentricities for Run2 in which planets have masses of 
nti = 10 M© (red line) and nto = 5 M© (blue line). The bottom 
panel displays the evolution of both the resonant angle i/^ = 
4Ao - 3Ai - 3oji (black) associated with the 4:3 resonance and 
the apsidal angle Aoj = ajf - ojo (green). 
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Fig. 3. This figure shows the evolution of the semimajor axes 
and eccentricities for Run3 in which planets have masses of 
nti = 20 M© (green line) and nto = 10 M© (red line). The 
bottom panel displays the evolution of both the resonant angle 
i// = 3Ao - 2Ai - IcDi (black) associated with the 3:2 resonance 
and the apsidal angle Ao; = coi - cOo (green). 



model the initial planetary mass distribution decreases out- 
ward, so that the innermost planet has mass of mi = 15 M©. 
As noted already, the planets are spaced in orbital radius as- 
suming a mutual separation of 5 RmH- Thus the other planets 
with masses of = 12.5, 10, 7.5, 5 M© are located initially at 
Up = 2.1, 2.4, 2.8, 3.1 respectively. 

Fig. [9] shows the evolution of the semimajor axes and eccentric- 
ities of planets for this model. At the beginning of the simula- 
tion, all the planets migrate inward as a consequence of type I 
migration, with a migration rate decreasing as one moves from 
the innermost planet to the outermost one. At ^ ~ 7 x 10^, the 
innermost and most rapidly migrating 15 M© planet reaches 
the edge of the inner cavity located at r ~ 1 .2. As expected, this 
body remains trapped at this location until the inwardly migrat- 
ing 12.5 M© planet approaches and they enter the 4:3 resonance 
at ^ ~ 1.1 X 10^. From this time the evolution of these two plan- 
ets is similar to that of pairs of planets with q > I described 
in Section [T2l the planets reach a quasi equilibrium state such 
that they evolve on non migrating orbits with ai ~ 1.2 and 
a2 ~ 1.4, and with eccentricities remaining almost constant. 
This lasts until the 10 M© planet enters the 4:3 resonance with 
the second body at ^ ~ 1.5x10"^. Again, the third planet tends to 
push the innermost planets inward, but the corotation torques 
exerted on the 15 M© planet are able to counterbalance this 



eff'ect and the migration of this three-planet system is stalled. 
A similar process occurs each time a migrating planet is reso- 
nantly captured by the bodies located inside its orbit. Over time 
we find that some planets slip from one resonance to another, 
but the system remains globally stable during these episodes. 
For example, the fifth planet with mass ms = 5 M© slips from 
the 4:3 resonance with the 7.5 M© planet to the 5:4 resonance at 
^ ~ 4.5x 10"^ and finally enters the 6:5 resonance at ^ ~ 5.2x 10"^. 
The final outcome of the simulation is a system forming a se- 
ries of resonances between adjacent bodies with each of them 
evolving on a non migrating orbit. 

Fig [TOldisplays the resonant angles j/^i = (p-\-l)Ao-pAi-0Ji 
and i//2 = (p l)/^o - p^i - oJo corresponding to the (/? -h 1) : 
p commensurabilities that form between each pair of adja- 
cent bodies. We see that all commensurabilities that form are 
first order resonances, in agreement with results obtained by 
Cress well & Nelson (2006). This simulation suggests that in 
a circumbinary disc, corotation torques exerted at the edge of 
the inner cavity provide an efficient mechanism against type 
I migration for a swarm of planets, and that resonant capture 
prevents close encounters, scattering and collisions when the 
initial planetary mass distribution decreases as a function of 
orbital radius. 
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Fig. 10. This figure shows the resonances which are established between adjacent bodies at the end of the simulation correspond- 
ing to Model 1. Planets are labelled from 1 to 5, with 1 being the innermost planet and 5 being the outermost one. 
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Fig. 4. This figure shows the evolution of the semimajor axes 
and eccentricities for Run4 in which planets have masses of 
Mi = 20 Me (green line) and = 5 Me (blue line). 

4.2. Model2 

In this model, the initial configuration of the system is such that 
moving from the innermost planet to the outermost one, planets 
have masses ofm^ = 15, 7.5, 12.5, 5 and 10 Me respectively. 
A snapshost of the disc surface density at the beginning of the 
simulation is presented in the left panel of Fig. [121 and the evo- 



Fig. 5. The upper panel shows the evolution of the torques ex- 
erted by the disc on the planets for Run4. The green line corre- 
sponds to the torques exerted on the 20 Me planet and the blue 
line corresponds to the torques exerted on the 5 Me planet. The 
lower panel displays the azimutal average of the surface den- 
sity at the time shown above the plot as well as the position of 
the planets at the same time. 

lution of the semimajor axes and eccentricities of planets for 
this model is illustrated in Fig. [TT] Once again, the innermost 
15 Me planet rapidly drifts toward the edge of the inner cavity 
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Fig. 6. This figure shows the evolution of the semimajor axes 
and eccentricities for Run5 in which planets have masses of 
nti = 5 M® (blue line) and = 10 M® (red line). The bottom 
panel displays the evolution of both the resonant angle i// = 
6Ao - 5Ai - 5(jJi (black) associated with the 6:5 resonance and 
the apsidal angle /^co = oji - cOo (green). 
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Fig. 7. This figure shows the evolution of the semimajor axes 
and eccentricities for Run6 in which planets have masses of 
Mi = 10 M© (red line) and = 20 M© (green line). The 
bottom panel displays the evolution of both the resonant angle 
i// = 6Ao - 5Ai - Soji (black) associated with the 6:5 resonance 
and the apsidal angle Aoj = oji - cOo (green). 



where it remains trapped. In comparison with Model 1, adja- 
cent bodies have significantly diff'erent masses which leads to 
a stronger diff'ential migration between them. This causes the 
12.5 M© and 10 M© planets to rapidly catch up with the 7.5 
M© and 5 M© planets respectively, leading to the formation of 
resonances. Here, the 12.5 M© planet first enters the 6:5 res- 
onance with the 7.5 M© body and then slips into the 7:6 res- 
onance. Fig. [TT] shows that the excitation of eccentricities due 
to this resonant interaction, and the influence of surrounding 
planets, leads ultimately to a collision between these two bod- 
ies at ^ ~ 1.7 x 10^. This merger forms a 20 M© planet, which 
then migrates inward until it becomes locked stably into the 
5:4 resonance with the 15 M© planet. At ^ ~ 3 x 10"^ a new col- 
lision occurs between the 5 and 10 M© planets. Earlier, these 
two planets passed through a sequence of difl'erent resonances 
and were in the 7:6 resonance just before this collision occured. 
Again, this newly formed 15 M© planet migrates until it catches 
with the 20 M© body. 

The right panel of Fig. [12] shows a snapshot of the disc sur- 
face density at the end of the simulation. We see that the final 
state of the system consists of only three bodies evolving in 
the disc. Moving from the innermost planet to the outermost 
one, these have masses of 15, 20 and 15 M©. Once again, the 
corotation torques exerted at the cavity edge prevent inward 



migration, leaving each surviving body in resonance with its 
neighbours. Fig.[T3]displays the resonant angles corresponding 
to the commensurabilities that form between each pair. These 
are the 5:4 resonance for the first pair and the 4:3 resonance 
for the second one. This is in good agreement with the simula- 
tions performed by Cresswell & Nelson (2006) who studied the 
evolution of a swarm of planets embedded in a protoplanetary 
disc, and who found that the 4:3, 5:4, 6:5 and 7:6 resonances 
are most favoured. 

4.3. Models 

In this last model, we consider planets with masses nip = 5, 7.5, 
10, 12.5 and 15 M© located initially at ap = 1.8, 2, 2.3, 2.7 and 
3.1 respectively. 

The evolution of the semimajor axes and eccentricities of plan- 
ets for Model3 are shown in Fig. [141 At the beginning of the 
simulation, each planet migrates inward, but the initial configu- 
ration of this model is such that the orbits of two adjacent bod- 
ies rapidly converge, leading to the formation of resonances. 
At ^ ~ 3 X 10^, the system has already evolved to a state 
where each body is in resonance with its closest neighbours, 
which occurs here before the innermost planet reaches the in- 
ner cavity. Moving from the inner planet to the outer one, the 
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Fig. 12. This figure shows, for Model2, snapshots of the disc surface density at times shown above the plots. In this figure, planets 
are represented by white circles. 
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Fig. 13. This figure shows the resonances which are established between adjacent bodies at the end of the simulation correspond- 
ing to Model2. Planets are labelled from 1 to 5, with 1 being the innermost planet and 5 being the outermost one. 



commensurabilities that form are respectively the 7:6, 6:5, 6:5 
and 5:4 resonances. The resonant interaction between the 5 and 
7.5 M© planets increases their eccentricities up to ep ~ 0.08. 
Eccentricity growth within the swarm leads to crossing orbits, 
and subsequent interactions cause the breaking of resonances. 
This leads to a collison between the 5 and 7.5 M© planets at 
r ~ 1.8 X 10"^. At r ~ 2.4 X 10"^ a colHsion resulting in a 22.5 M© 
planet also occurs between the 10 and 12.5 M© bodies which 
were locked into the 7:6 resonance prior to this event. 

The final state of the system consists of three planets with 
masses of = 12.5, 22.5 and 15 M© respectively located at 
Gp = 1.1, 1.3 and 1.5. The upper panel of Fig. [15] shows that 
the 5:4 resonance is clearly established between the first pair of 
planets. At the end of the simulation, these two bodies evolve 
on non migrating orbits, whereas the outermost body migrates 
outward very slightly, indicating that the latter is not in reso- 
nance with its interior neighbour. The disc surface density pro- 
file at r ~ 4 X 10"^ is displayed in the bottom panel of Fig. \15\ 



It shows that the outermost planet is located in a region where 
the disc has a large positive surface density gradient, close to 
the outer edge of the partial gap formed by the 22.5 M© body. 
This results in large positive corotation torques, leading to the 
observed outward migration of the outermost planet. This pro- 
cess is likely to operate until the planet reaches a fixed point 
located near the gap edge, where the total torque (corotation 
plus Lindblad) cancels (Masset et al. 2006). This result is con- 
sistent with that found in Section [T2l where the inner and outer 
planet masses were 20 and 5 M©, respectively. 

4.4. Long term evolution 

We now turn to the question of the long-term evolution of 
the planetary systems obtained in the five-planet simulations. 
Because the interaction with the gas disc tends to damp ec- 
centricities, it is necessary to examine the dynamical stability 
of the planets after the disc dispersal to establish long term 
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Fig. 8. This figure shows the evolution of the semimajor axes 
and eccentricities for Run6 in which planets have masses of 
rrii = 5 Me (blue line) and = 20 Me (green line). 
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Fig. 9. This figure shows the evolution of the semimajor axes 
and eccentricities of planets for Model 1. Moving from the in- 
nermost planet to the outermost one, planets have masses of 15 
(black line), 12.5 (red line), 10 (green line), 7.5 (blue line) and 
5 Me (yellow line). 



stability. Each of the previous simulations was restarted at a 
point corresponding to the end of the run, but with the gas 
surface density decaying exponentially with an e-folding time 
tdec = 2 X 10^. Once these systems had evolved for ~ 10^ bi- 
nary orbits, by which time the surface density in the discs had 
decreased by a factor of ~ 10^, we continued the simulations 
with a pure N-body code, ignoring any residual eff'ects of the 
remaining gas. For each of the five-planet models the results of 
this procedure are presented in Fig.O which displays the time 
evolution of the orbital radii of the planets. 

In Model 1, the eccentricity growth resulting from the disc 
dispersal gives rise, at the beginning of the simulation, to nu- 



Fig. 11. This figure shows the evolution of the semimajor axes 
and eccentricities of planets for Model2. Moving from the in- 
nermost planet to the outermost one, planets have masses of 15 
(black line), 7.5 (red line), 12.5 (green line), 5 (blue line) and 
10 Me (yellow line). 
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Fig. 14. This figure shows the evolution of the semimajor axes 
and eccentricities of planets for Model3. Moving from the in- 
nermost planet to the outermost one, planets have masses of 5 
(black line), 7.5 (red line), 10 (green line), 12.5 (blue line) and 
15 Me (yellow line). 



merous scattering events that eventually lead to collisions. At 
time t - 8 X 10^ a system of two planets with masses of 27.5 
and 22.5 Me remains, but these merge at ^ ~ 1.1x10^. The final 
state of the system is then a 50 Me planet orbiting at r ~ 1.5. 

A similar outcome is obtained in Model3 in which the long- 
term evolution resulted in a 15 Me planet evolving in a high- 
eccentricity orbit with a semi-major axis of Gp ~ 2. At earlier 
times, the increase in eccentricities following disc dispersal led 
to a collision between the 12.5 and 22.5 Me bodies, thereby 
forming a new 35 Me planet. At t ~ 2 x 10^, the latter is 
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Fig. 15. The upper panel shows, for Model3, the resonant an- 
gles (black) and i/^2 (green) associated with the resonance 
that forms between the two innermost planets. The bottom 
panel displays the disc surface density profile at the time shown 
above the plot. The semimajor axes of the planets at that time 
are represented by black circles. 

observed to undergo a close encounter with the cental binary, 
leading to this body being completely ejected from the system. 
Interestingly, the three-planet system in Model2 appears to 
be dynamically stable over long time scales, with the plan- 
ets maintaining their commensurabilities. This indicates that 
multiplanet resonant systems could potentially be found in cir- 
cumbinary discs, where the existence of the resonance helps to 
maintain the stability of the system. 

5. Summary and conclusion 

In this paper we have presented the results of hydrodynamic 
simulations aimed at studying the evolution of multiple planets 
embedded in a circumbinary disc. 

We first focused on a system consisting of a pair of planets in- 
teracting with each other. We assumed that one body is trapped 
at the edge of the inner cavity formed by the binary, while the 
other migrates inward from outside the orbit of the innermost 
body. Our calculations show diff'erent outcomes, depending on 
the planet mass ratio q = mi/mo. 

i) For models with q = 1 , the simulations indicate that the sys- 
tem evolves toward a quasi equilibrium state with both plan- 
ets trapped in a mean motion resonance and evolving on non 
migrating orbits. This occurs because the positive corotation 
torques exerted on the innermost planet counterbalance the 
negative corotation torques exerted on the outermost one. 

ii) For models with q > I, most of the runs resulted in a similar 
mode of evolution. Occasionally however, the final fate of the 
sytem was such that the two bodies are in close vicinity to one 
another, have stopped migrating, but are not in resonance. This 
arises when the mass of the innermost planet is high enough 
to open a partial gap in the disc, such that the migration of the 
outer planet is stopped at the edge of the large cavity formed 
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Fig. 16. This figure shows time evolution of the orbital radii of 
the planets for Model 1, Model2 and Model3 as a consequence 
of the disc dispersal. It corresponds to the results of hydro plus 
N-body simulations (see text for details). 



by both the binary and the innermost planet, 
iii) For most of the models with q < 1 , the planets involved in 
the simulations underwent diff'erent dynamical processes such 
as scattering or orbital exchange. When orbital exchange oc- 
curs, the final state of the system is a stable mean motion res- 
onance with the more massive planet now being the innermost 
one. Scattering and orbital exchange occurs once a first order 
resonance is established between the planets. This drives up 
the eccentricity of the inner planet leading to a close encounter 
with either the outer planet or the binary. A close encounter 
with the binary leads to ejection. 

Having examined the two-planet problem, we then focused 
on more "realistic" systems composed of five planets with 
masses of 5, 7.5, 10, 12.5 and 15 M© embedded in a circumbi- 
nary disc. We performed three simulations, with the planets be- 
ing placed in diff'erent initial orbital configurations. 

In general terms, the evolution of such a system proceeds 
as follows. Each planet migrates inward until it is captured into 
resonance with an adjacent body. This occurs either because 
the innermost planet has reached the edge of the inner cav- 
ity, where it remains trapped, or because two adjacent bodies 
migrate diff'erentially. In some simulations, resonant interac- 
tion and the associated eccentricity growth resulted in close 
encounters and collisions between bodies, resulting in the for- 
mation of a more massive planet. At later times, once the more 
chaotic phase of evolution is over, the system behaves more 
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quiescently, and generally reaches a stationary state with each 
remaining body in resonance with its closest neighbours on non 
migrating orbits due to planet trapping at the cavity edge. In 
two of the three simulations, we found that disc dispersal led 
to strong scattering and collisions between bodies, resulting in 
eventually only one planet being left in the system. In the re- 
maining case, we found that three planets could exist stably in 
a three body mean motion resonance over very long times. 

The results of our simulations indicate that resonant sys- 
tems of planets with masses in the ~ few-Earth mass range 
may be common in close binary systems. Indeed, the inner edge 
of a circumbinary disc plays the role of a barrier, which raises 
the possibility that two bodies become resonantly trapped. This 
resonance trapping, however, usually prevents the two planets 
from approaching each other so that they are unable to merge 
to form a larger body. When more than two bodies are intro- 
duced we found that the resulting scattering and close encoun- 
ters could lead to planetary growth, however, but we have been 
unable to run sufficient numbers of simulations to provide reli- 
able statistical information about the range of outcomes. 

A number of outstanding issues remain when considering 
planet formation in circumbinary discs. The first is understand- 
ing where the planetary cores can form due to planetesimal ac- 
cretion. The formation of an eccentric disc around a binary sys- 
tem, and the presence of the binary itself, can increase the ve- 
locity dispersion to a value where accretion does not occur ex- 
cept far out in the disc. The second is understanding the growth 
of cores once they have formed, either through continued plan- 
etesimal accretion, or through giant impacts. The migration of 
the core toward the inner cavity may actually place it in a re- 
gion where planetesimal accretion is reduced due to increased 
velocity dispersion, and at present there have been no simula- 
tions of giant impact growth which include the gas disc. If a 
giant core can be formed, however, then it is likely that a gas 
giant planet will result. As it grows, a giant planet should leave 
the edge of the cavity and undergo type II migration (e.g. Lin & 
Papaloizou 1993; Nelson & al. 2000), until it becomes trapped 
into the 4:1 resonance with the binary. This appears to be the 
most likely fate of a giant planet migrating in a circumbinary 
disc (Nelson 2003), although trapping is not necessarily perma- 
nent and simulations indicate a finite probablity of the planet 
being ejected from the system due to close encounters with the 
binary. The evolution of the planet during its growth from a 
core into a gas giant, however, has not been explored, and this 
will be the subject of a future paper. 
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